Navigation system with apparatus for detecting accuracy failures

ABSTRACT

A navigation system for a vehicle having a receiver operable to receive a plurality of signals from a plurality of transmitters includes a processor and a memory device. The memory device has stored thereon machine-readable instructions that, when executed by the processor, enable the processor to determine a set of error estimates corresponding to pseudo-range measurements derived from the plurality of signals, determine an error covariance matrix for a main navigation solution using ionospheric-delay data, and, using a parity space technique, determine at least one protection level value based on the error covariance matrix.

CROSS REFERENCE TO RELATED APPLICATION

This application claims priority from U.S. Provisional Appl. No.61/012,303 entitled “RAIM WITH SPATIALLY CORRELATED IONOSPHERIC ERRORS”filed Dec. 7, 2007, which is incorporated by reference as if fully setforth herein.

BACKGROUND OF THE INVENTION

Conventional RAIM algorithms may be based on either a weighted orun-weighted least squares solution where the errors in each satellite'spseudo-range measurement are uncorrelated with the errors in the othersatellites' pseudo-range measurements.

However, the ionospheric error (which can be the dominant error source)in each satellite's pseudo-range is, in fact, highly correlated withthat of each of the other satellites. By ignoring this correlation, thecomputed Horizontal Protection Limit (HPL) which bounds the horizontalposition error is much larger than necessary. As a result theavailability of GPS to do a low Required Navigation Performance (RNP)approach suffers.

SUMMARY OF THE INVENTION

In an embodiment of the invention, a navigation system for a vehiclehaving a receiver operable to receive a plurality of signals from aplurality of transmitters includes a processor and a memory device. Thememory device has stored thereon machine-readable instructions that,when executed by the processor, enable the processor to determine a setof error estimates corresponding to pseudo-range measurements derivedfrom the plurality of signals, determine an error covariance matrix fora main navigation solution using ionospheric-delay data, and, using aparity space technique, determine at least one protection level valuebased on the error covariance matrix.

BRIEF DESCRIPTION OF THE DRAWINGS

Preferred and alternative embodiments of the present invention aredescribed in detail below with reference to the following drawings.

FIG. 1 shows a navigation system incorporating embodiments of thepresent invention; and

FIG. 2 shows a process according to an embodiment of the invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

FIG. 1 shows a radio navigation system incorporating the teachings of anembodiment of the present invention. The system includes severaltransmitters 1-N and user set 12. Transmitters 1-N, in the preferredembodiment, are a subset of the NAVSTAR GPS constellation of satellitetransmitters, with each transmitter visible from the antenna of user set12. Transmitters 1-N broadcast N respective signals indicatingrespective transmitter positions and signal transmission times to userset 12.

User set 12, mounted to an aircraft (not shown), includes receiver 14,processor 16, and a memory device, such as processor memory 18. Receiver14, preferably NAVSTAR GPS compatible, receives the signals, extractsthe position and time data, and provides pseudorange measurements toprocessor 16. From the pseudorange measurements, processor 16 derives aposition solution for the user set 12. Although the satellites transmittheir positions in World Geodetic System of 1984 (WGS-84) coordinates, aCartesian earth-centered earth-fixed system, the preferred embodimentdetermines the position solution in a local reference frame L, which islevel with the north-east coordinate plane and tangential to the Earth.This frame choice, however, is not critical, since it is well-understoodhow to transform coordinates from one frame to another.

Processor 16 also uses the pseudorange measurements to detect satellitetransmitter failures and to determine a worst-case error, or protectionlimit, both of which it outputs with the position solution to flightmanagement system 20. Flight management system 20 compares theprotection limit to an alarm limit corresponding to a particularaircraft flight phase. For example, during a pre-landing flight phase,such as nonprecision approach, the alarm limit (or allowable radialerror) is 0.3 nautical miles, but during a less-demanding oceanic flightphase, the alarm limit is 2-10 nautical miles. (For more details onthese limits, see RTCA publication DO-208, which is incorporated hereinby reference.) If the protection limit exceeds the alarm limit, theflight management system, or its equivalent, announces or signals anintegrity failure to a navigational display in the cockpit of theaircraft. The processor also signals whether it has detected anysatellite transmitter failures.

An embodiment of the invention models the correlation of the ionosphericerrors between each pair of satellites as a function of the distancebetween their ionospheric pierce points. The closer the pierce points,the higher the correlation. The root-mean-square (RMS) uncertainty (orsigma) of each satellite's pseudo-range measurement is computed usingthe ionospheric variance model defined in DO-229D, Appendix J. Using thecomputed correlation coefficients and the sigma for each satellite, theionospheric measurement error covariance matrix is formed. The remainingerrors (satellite clock and ephemeris, tropospheric, multi-path andreceiver noise) are assumed to be uncorrelated. Thus, the combinedmeasurement error covariance matrix for these error sources is diagonal.These two matrices are added to form the total measurement errorcovariance matrix. This matrix is then inverted to form the weightingmatrix for the least squares solution. Fault detection and exclusion canthen be performed and the various protection levels such as thehorizontal protection level (HPL), vertical protection level (VPL),horizontal exclusion level (HEL), and vertical exclusion level (VEL)computed based on the methods of solution separation previouslydescribed in U.S. Pat. Nos. 5,760,737 and 6,639,549, each of which areincorporated by reference as if fully set forth herein.

FIG. 2 illustrates a process 200, according to an embodiment of theinvention, that can be implemented in the radio navigation systemillustrated in FIG. 1. The process 200 is illustrated as a set ofoperations or steps shown as discrete blocks. The process 200 may beimplemented in any suitable hardware, software, firmware, or combinationthereof. As such the process 200 may be implemented incomputer-executable instructions that can be transferred from oneelectronic device to a second electronic device via a communicationsmedium. The order in which the operations are described is not to benecessarily construed as a limitation.

Referring to FIG. 2, at a step 210, the processor 16 computes the sigma(error) values on pseudo-range and measurements.

At a step 220, the processor 16 determines the measurement matrix. Thetrue vector of pseudo-range residuals Δρ is related to the incrementalposition/time solution vector Δx (distance from the positionlinearization point) as follows:

Δρ=ρ−{circumflex over (ρ)}=HΔx  (1)

where H is the measurement matrix and is given by:

$H = {{\begin{bmatrix}{- {LOS}_{1\; x}} & {- {LOS}_{1y}} & {- {LOS}_{1\; z}} & 1 \\{- {LOS}_{2\; x}} & {- {LOS}_{2\; y}} & {- {LOS}_{2\; z}} & 1 \\\vdots & \vdots & \vdots & \vdots \\{- {LOS}_{N\; x}} & {- {LOS}_{Ny}} & {- {LOS}_{Nz}} & 1\end{bmatrix}\begin{bmatrix}{LOS}_{ix} \\{LOS}_{iy} \\{LOS}_{iz}\end{bmatrix}} = \begin{matrix}{{{Line}\text{-}{of}\text{-}{sight}\mspace{14mu} {unit}\mspace{14mu} {vector}\mspace{14mu} {pointing}}\mspace{14mu}} \\{{from}\mspace{14mu} {the}\mspace{14mu} {user}\mspace{14mu} {to}\mspace{14mu} {satellite}\mspace{14mu} i}\end{matrix}}$Δ x = x − x̂ = true  position/clock  bias-linearization  point 

At a step 230, the processor 16 computes the Error Covariance Matrix.The vector of measured pseudo-range residuals Δ{tilde over (ρ)} is thetrue pseudo-range residual vector plus the vector of residual errors δρand is thus:

$\begin{matrix}\begin{matrix}{{\Delta \overset{\sim}{\rho}} = {\rho + {\delta \; \rho} - \hat{\rho}}} \\{= {\overset{\sim}{\rho} - \hat{\rho}}} \\{= {{H\; \Delta \; x} + {\delta \; \rho}}}\end{matrix} & (2)\end{matrix}$

The processor 16 designates the post-update estimate of Δx asΔ{circumflex over (x)}. Then, the processor 16 can define the vector ofpost-update measurement residuals as:

ξ=Δ{tilde over (ρ)}−HΔ{circumflex over (x)}  (3)

Each post-update measurement residual is the difference between themeasured pseudo-range residual and the predicted pseudo-range residualbased on the post-update estimate Δ{circumflex over (x)}.

A “weighted least-squares solution” can be determined by the processor16 by finding the value of Δ{circumflex over (x)} which minimizes theweighted sum of squared residuals. Thus, the processor 16 may minimize:

ξ^(T) Wξ=(Δ{tilde over (ρ)}−HΔ{circumflex over (x)})^(T) W(Δ{tilde over(ρ)}−HΔ{circumflex over (x)})  (4)

where W is an appropriate weighting matrix. The weighting matrixgenerally chosen is one which normalizes the residuals based on theuncertainty of each pseudo-range measurement. Thus, the processor 16yields:

$\begin{matrix}{{W = \begin{bmatrix}{1/\sigma_{1}^{2}} & 0 & \ldots & 0 \\0 & {1/\sigma_{2}^{2}} & \ldots & 0 \\\vdots & \vdots & \ddots & \vdots \\0 & 0 & \ldots & {1/\sigma_{N}^{2}}\end{bmatrix}},\mspace{14mu} \begin{matrix}{{{assuming}\mspace{14mu} {uncorrelated}}\mspace{14mu}} \\{measurements}\end{matrix}} & (5)\end{matrix}$

which represents the inverse of the pseudo-range measurement errorcovariance matrix assuming each pseudo-range error is uncorrelated withthe others.

However, the vertical ionospheric delay component of each pseudo-rangeerror is highly correlated with the others. If this correlation isknown, then the processor 16 can take advantage of that knowledge byusing the true pseudo-range measurement error covariance matrix R. Theweighting matrix then becomes

$\begin{matrix}{W = {R^{- 1} = \begin{bmatrix}\sigma_{1}^{2} & {E\left\lbrack {\delta \; \rho_{1}\delta \; \rho_{2}} \right\rbrack} & \ldots & {E\left\lbrack {\delta \; \rho_{1}\delta \; \rho_{N}} \right\rbrack} \\{E\left\lbrack {\delta \; \rho_{1}\delta \; \rho_{2}} \right\rbrack} & \sigma_{2}^{2} & \ldots & {E\left\lbrack {\delta \; \rho_{2}\delta \; \rho_{N}} \right\rbrack} \\\vdots & \vdots & \ddots & \vdots \\{E\left\lbrack {\delta \; \rho_{1}\delta \; \rho_{N}} \right\rbrack} & {E\left\lbrack {\delta \; \rho_{2}\delta \; \rho_{N}} \right\rbrack} & \ldots & \sigma_{N}^{2}\end{bmatrix}^{- 1}}} & (6)\end{matrix}$

The value of Δ{circumflex over (x)} that minimizes (4) is determined bytaking the derivative, setting it equal to zero, and solving forΔ{circumflex over (x)}. This yields:

$\begin{matrix}\begin{matrix}{{{\Delta \hat{x}} = {\left( {H^{T}{WH}} \right)^{- 1}H^{T}W\; \Delta \; \overset{\sim}{\rho}}}\;} \\{= {S\; \Delta \; \overset{\sim}{\rho}}}\end{matrix} & (7)\end{matrix}$

where the processor 16 has defined the weighted least-squares solutionmatrix S as:

S=(H^(T)WH)⁻¹H^(T)W  (8)

Altitude Aiding

Barometric altitude can be used by the processor 16 to augment the GPSpseudo-range measurements. If it is used, the measurement matrix is thenaugmented as follows

$\begin{matrix}{H = \begin{bmatrix}{- {LOS}_{1\; x}} & {- {LOS}_{1y}} & {- {LOS}_{1\; z}} & 1 \\{- {LOS}_{2\; x}} & {- {LOS}_{2\; y}} & {- {LOS}_{2\; z}} & 1 \\\vdots & \vdots & \vdots & \vdots \\\begin{matrix}{- {LOS}_{N\; x}} \\0\end{matrix} & \begin{matrix}{- {LOS}_{Ny}} \\0\end{matrix} & \begin{matrix}{- {LOS}_{Nz}} \\{- 1}\end{matrix} & \begin{matrix}1 \\0\end{matrix}\end{bmatrix}} & (9)\end{matrix}$

This measurement matrix assumes that the incremental position vector(the first 3 elements) within Δx are given in local-level coordinates(with the z axis down). The line-of-sight (LOS) elements then must alsobe expressed in the local-level coordinates. The weighting matrix isalso augmented as follows

$\begin{matrix}{W = {R^{- 1} = \begin{bmatrix}\sigma_{1}^{2} & {E\left\lbrack {\delta \; \rho_{1}\delta \; \rho_{2}} \right\rbrack} & \ldots & {E\left\lbrack {\delta \; \rho_{1}\delta \; \rho_{N}} \right\rbrack} & 0 \\{E\left\lbrack {\delta \; \rho_{1}\delta \; \rho_{2}} \right\rbrack} & \sigma_{2}^{2} & \ldots & {E\left\lbrack {\delta \; \rho_{2}\delta \; \rho_{N}} \right\rbrack} & 0 \\\vdots & \vdots & \ddots & \vdots & \vdots \\{E\left\lbrack {\delta \; \rho_{1}\delta \; \rho_{N}} \right\rbrack} & {E\left\lbrack {\delta \; \rho_{2}\delta \; \rho_{N}} \right\rbrack} & \ldots & \sigma_{N}^{2} & 0 \\0 & 0 & \ldots & 0 & \sigma_{baro}^{2}\end{bmatrix}^{- 1}}} & (10)\end{matrix}$

Computing the Measurement Covariance Matrix

There are multiple methods that can be employed to determine themeasurement error covariance matrix. In the case of a Kalman filterapplication, the temporal behavior (time-correlation) of the ionosphericdelays may be modeled. The spatially correlated ionospheric error forsatellite i can be modeled as a weighted sum of three independentnormalized (sigma=1.0) Gaussian random errors scaled by the nominal ionosigma value for that satellite as follows:

δρ_(iono) _(—) _(i)=σ_(iono) _(—) _(i) k _(iono) _(—) _(i) ^(T) x_(ref)  (11)

where x_(ref) is a 3×1 vector of independent Gaussian random errors withzero mean and a variance of 1. The weighting vector k_(iono) _(—) _(i)^(T) is determined by the processor 16 by first defining three gridpoints on the thin shell model of the ionosphere (at a height of 350 km)equally spaced in azimuth at a great circle distance of 1500 km from theuser. The processor 16 may then define a 3×1 vector x_(grid) ofnormalized delays at these points. The delays at gridpoints i and j maybe spatially correlated with each other based on the great circledistances between them according to:

E[x_(grid) _(—) _(i)x_(grid) _(—) _(i)^(T)]=1−(1−e^(−d)grid/diono)²  (12)

where

d_(grid) _(—) _(i,grid) _(—) _(j)=great circle distance between gridpoint i and grid point j

d_(iono)=correlation distance of the ionospheric delay=4000 km  (13)

Using that relationship the processor 16 can form a 3×3 covariancematrix P_(grid) which describes the correlations between each of thegrid points:

P_(grid)=E[x_(grid)x_(grid) ^(T)]

If the delay processes that exist at these grid points are a certainlinear combination of the reference independent Gaussian random errors,then they will have the desired spatial and temporal correlation. Theprocessor 16 may assume that the desired linear combination is obtainedby using a 3×3 upper-triangular mapping matrix U_(grid) as follows:

x_(grid)=U_(grid)x_(ref)  (14)

The grid covariance matrix is then:

$\begin{matrix}\begin{matrix}{P_{grid} = {{E\left\lbrack {x_{grid}x_{grid}^{T}} \right\rbrack} = {U_{grid}{E\left\lbrack {x_{ref}x_{ref}^{T}} \right\rbrack}U_{grid}^{T}}}} \\{= {U_{grid}U_{grid}^{T}}}\end{matrix} & (15)\end{matrix}$

Therefore, the mapping matrix U_(grid) can be formed by the processor 16simply by factoring the covariance matrix P_(grid). Since the geometryof the three gridpoints is fixed, the covariance matrix P_(grid) isconstant and can thus be pre-computed by the processor 16. Now theprocessor 16 can choose a linear combination of the three grid-pointdelays that yields a normalized delay at the pierce-point of thesatellite i such that the proper spatial correlation with the three gridpoints (and thus, presumably, each of the other satellites) is achievedas follows:

δρ_(norm) _(—) _(iono) _(—) _(i) =k _(sat) _(—) _(i) _(—) _(grid) ^(T) x_(grid)  (16)

where

k_(sat) _(—) _(i) _(—) _(grid)=3-vector of weighting factors

δρ_(norm) _(—) _(iono) _(—) _(i)=normalized delay at the satellitepierce-point

The satellite pseudo-range delay may be correlated to the delay at thek^(th) grid point according to:

where

d_(sat) _(—) _(i,grid) _(—) _(k)=great circle distance between thesatellite pierce point and the grid point  (17)

d_(iono)=correlation distance of the nominal ionospheric delay  (18)

The 1×3 covariance matrix P_(sat) _(—) _(i) _(—) _(grid), which definesthe correlations between satellite i and each of the grid points, is

$\begin{matrix}\begin{matrix}{P_{{sat\_ i}{\_ grid}} = {{E\left\lbrack {\delta \; \rho_{norm\_ i}x_{grid}^{T}} \right\rbrack} = {E\left\lbrack {K_{{sat\_ i}{\_ grid}}^{T}x_{grid}x_{grid}^{T}} \right\rbrack}}} \\{= {K_{{sat\_ i}{\_ grid}}^{T}P_{grid}}}\end{matrix} & (19)\end{matrix}$

Therefore the weighting vector k_(sat) _(—) _(grid) can be found by theprocessor 16 as follows:

k_(sat) _(—) _(i) _(—) _(grid) ^(T)=P_(sat) _(—) _(i) _(—)_(grid)P_(grid) ⁻¹  (20)

Combining (14) and (16), the processor 16 can obtain the normalizedvertical delays directly from the three independent reference delays asfollows:

$\begin{matrix}\begin{matrix}{{\delta \; \rho_{{norm\_ iono}{\_ i}}} = {K_{{sat\_ i}{\_ grid}}^{T}U_{grid}x_{ref}}} \\{= {K_{iono\_ i}^{T}x_{ref}}}\end{matrix} & (21)\end{matrix}$

Thus, the weighting vector is:

k_(iono) _(—) _(i) ^(T)=k_(sat) _(—) _(i) _(—) _(grid)U_(grid)  (22)

The processor 16 can form a vector of N normalized pseudo-range ionodelays from (21) as follows:

$\begin{matrix}\begin{matrix}{{\delta \; \rho_{{norm\_ iono}{\_ i}}} = {\begin{bmatrix}K_{{sat\_}1{\_ grid}}^{T} \\K_{{sat\_}2{\_ grid}}^{T} \\\vdots \\K_{{{sat\_}N}{\_ grid}}^{T}\end{bmatrix}U_{grid}x_{ref}}} \\{= {K_{sat\_ grid}U_{grid}x_{ref}}}\end{matrix} & (23)\end{matrix}$

The actual (non-normalized) delay along the line of sight can beobtained by the processor 16 by scaling the normalized delay by thesigma value for that satellite based on the geomagnetic latitude of thepierce-point and obliquity factor as defined in DO-229. In vector form,the processor 16 yields:

$\begin{matrix}\begin{matrix}{{\delta \; \rho_{iono}} = {\begin{bmatrix}\sigma_{{iono\_}1} & 0 & \ldots & 0 \\0 & \sigma_{{iono\_}2} & \ddots & \vdots \\\vdots & \ddots & \ddots & 0 \\0 & \ldots & 0 & \sigma_{iono\_ N}\end{bmatrix}\delta \; \rho_{norm\_ iono}}} \\{= {\begin{bmatrix}\sigma_{{iono\_}1} & 0 & \ldots & 0 \\0 & \sigma_{{iono\_}2} & \ddots & \vdots \\\vdots & \ddots & \ddots & 0 \\0 & \ldots & 0 & \sigma_{iono\_ N}\end{bmatrix}K_{sat\_ grid}U_{grid}x_{ref}}} \\{= {\Gamma \; K_{sat\_ grid}U_{grid}x_{ref}}}\end{matrix} & (24)\end{matrix}$

The ionospheric delay error covariance matrix may be defined as:

$\begin{matrix}\begin{matrix}{R_{iono} = {E\left\lbrack {{\delta\rho}_{iono}{\delta\rho}_{iono}^{T}} \right\rbrack}} \\{= {\Gamma \; K_{sat\_ grid}U_{grid}{E\left\lbrack {x_{ref}x_{ref}^{T}} \right\rbrack}U_{grid}^{T}\; K_{sat\_ grid}^{T}\Gamma^{T}}} \\{= {\Gamma \; K_{sat\_ grid}P_{grid}K_{sat\_ grid}^{T}\Gamma^{T}}}\end{matrix} & (25) \\{R_{iono} = \left\lbrack \begin{matrix}\sigma_{{iono\_}\; 1}^{2} & {E\left\lbrack {{\delta\rho}_{{iono\_}\; 1}{\delta\rho}_{{iono\_}\; 2}} \right\rbrack} & \ldots & {E\left\lbrack {{\delta\rho}_{{iono\_}\; 1}{\delta\rho}_{{iono\_}\; N}} \right\rbrack} \\{E\left\lbrack {{\delta\rho}_{{iono\_}\; 1}{\delta\rho}_{{iono\_}\; 2}} \right\rbrack} & \sigma_{{iono\_}\; 2}^{2} & \cdots & {E\left\lbrack {{\delta\rho}_{{iono\_}\; 2}{\delta\rho}_{{iono\_}\; N}} \right\rbrack} \\\vdots & \vdots & \ddots & \vdots \\{E\left\lbrack {{\delta\rho}_{{iono\_}\; 1}{\delta\rho}_{{iono\_}\; N}} \right\rbrack} & {E\left\lbrack {{\delta\rho}_{{iono\_}2}{\delta\rho}_{{iono\_}\; N}} \right\rbrack} & \cdots & \sigma_{{iono\_}\; N}^{2}\end{matrix} \right\rbrack} & \;\end{matrix}$

The rest of the pseudo-range measurement errors are assumed to beuncorrelated with a composite one-sigma value denoted by σ_(other) _(—)_(i) for satellite i. For simplicity, the processor 16 can assume thatthe one-sigma value for each satellite is a constant six meters. Thetotal measurement error covariance matrix is then:

$\begin{matrix}{R = {W^{- 1} = {R_{iono} + \begin{bmatrix}\sigma_{{other\_}\; 1}^{2} & 0 & \ldots & 0 \\0 & \sigma_{{other\_}\; 2}^{2} & \ddots & \vdots \\\vdots & \ddots & \ddots & 0 \\0 & \ldots & 0 & \sigma_{{other\_}\; N}^{2}\end{bmatrix}}}} & (26)\end{matrix}$

In a snapshot RAIM approach, the correlations between satellites arecomputed directly without the use of a grid. Computing the correlationsbetween satellites directly may be both simpler and slightly moreaccurate.

Specifically, ionospheric error covariance may be modeled as a functionof the great circle distance between the pierce points along theionospheric shell (350 km above the earth's surface):

E[δρ _(iono) _(—) _(i)δρ_(iono) _(—) _(j)]=σ_(iono) _(—) _(i)σ_(iono)_(—) _(j)[1−(1−e ^(−d) ^(ij) ^(/d) ^(iono) )²]

where:

d_(ij)=great circle distance between pierce points for sats i and j

d_(iono)=de-correlation distance=4000 km

Ionospheric errors are highly correlated. As such:

W = R⁻¹ = (R_(iono) + R_(other))⁻¹ where:${R_{iono} = {{\left\lbrack \begin{matrix}\sigma_{{iono\_}\; 1}^{2} & {E\left\lbrack {{\delta\rho}_{{iono\_}\; 1}{\delta\rho}_{{iono\_}\; 2}} \right\rbrack} & \ldots & {E\left\lbrack {{\delta\rho}_{{iono\_}\; 1}{\delta\rho}_{{iono\_}\; N}} \right\rbrack} \\{E\left\lbrack {{\delta\rho}_{{iono\_}\; 1}{\delta\rho}_{{iono\_}\; 2}} \right\rbrack} & \sigma_{{iono\_}\; 2}^{2} & \cdots & {E\left\lbrack {{\delta\rho}_{{iono\_}\; 2}{\delta\rho}_{{iono\_}\; N}} \right\rbrack} \\\vdots & \vdots & \ddots & \vdots \\{E\left\lbrack {{\delta\rho}_{{iono\_}\; 1}{\delta\rho}_{{iono\_}\; N}} \right\rbrack} & {E\left\lbrack {{\delta\rho}_{{iono\_}2}{\delta\rho}_{{iono\_}\; N}} \right\rbrack} & \cdots & \sigma_{{iono\_}\; N}^{2}\end{matrix} \right\rbrack R_{other}} = \begin{bmatrix}\sigma_{{other\_}\; 1}^{2} & 0 & \ldots & 0 \\0 & \sigma_{{other\_}\; 2}^{2} & \ldots & 0 \\\vdots & \vdots & \ddots & \vdots \\0 & 0 & \ldots & \sigma_{{other\_}\; N}^{2}\end{bmatrix}}}{\; \mspace{290mu}}$

Error Covariance for the Weighted Least Squares Solution

At a step 240, the processor 16 computes a weighted least-squaressolution. The error in the post-updated solution is:

$\begin{matrix}\begin{matrix}{{\delta \hat{x}} = {{\Delta \hat{x}} - {\Delta \; x}}} \\{= {{\left( {H^{T}{WH}} \right)^{- 1}H^{T}W\; \Delta \overset{\sim}{\rho}} - {\Delta \; x}}}\end{matrix} & (27)\end{matrix}$

Substituting (2) into (27) yields:

$\begin{matrix}\begin{matrix}{{\delta \hat{x}} = {{\left( {H^{T}{WH}} \right)^{- 1}H^{T}{W\left( {{H\; \Delta \; x} + {\delta\rho}} \right)}} - {\Delta \; x}}} \\{= {{\Delta \; x} - {\left( {H^{T}{WH}} \right)^{- 1}H^{T}W\; {\delta\rho}} - {\Delta \; x}}} \\{= {\left( {H^{T}{WH}} \right)^{- 1}H^{T}W\; {\delta\rho}}} \\{= {S\; {\delta\rho}}}\end{matrix} & (28)\end{matrix}$

Thus, the solution matrix S maps the pseudo-range errors into thepost-updated solution error vector. The solution error covariance matrixmay be defined as:

$\begin{matrix}\begin{matrix}{P = {{E\left\lbrack {\delta \hat{x}\delta {\hat{x}}^{T}} \right\rbrack} - {{{SE}\left\lbrack {\delta\rho\delta\rho}^{T} \right\rbrack}S^{T}}}} \\{= {{SW}^{- 1}S^{T}}} \\{= {\left( {H^{T}{WH}} \right)^{- 1}H^{T}{WW}^{- 1}{{WH}\left( {H^{T}{WH}} \right)}^{- 1}}} \\{= \left( {H^{T}{WH}} \right)^{- 1}}\end{matrix} & (29)\end{matrix}$

The x and y horizontal position errors are statistically described bythe upper 2×2 portion of P. The major and minor axes of the horizontalposition error ellipse are equal to the square roots of the maximum andminimum eigenvalues of this 2×2 matrix and represent the one-sigmaerrors in the corresponding directions. Thus, the one-sigma error in theworst-case direction is given by:

$\begin{matrix}\begin{matrix}{\sigma_{horz\_ max} = \sqrt{\lambda_{\max}^{P{({{1\text{:}2},{1\text{:}2}})}}}} \\{= \sqrt{\frac{p_{11} + p_{22}}{2} + \sqrt{\left( \frac{p_{11} - p_{22}}{2} \right)^{2} + \left( p_{12} \right)^{2}}}}\end{matrix} & (30)\end{matrix}$

The one-sigma error in the vertical position is given by:

σ_(vert)√{square root over (p₃₃)}  (31)

Horizontal Figure of Merit is a conservative 95% fault-free error boundand may be computed by the processor 16 as the 2D RMS error from theerror covariance matrix

HFOM=2√{square root over (P(1,1)+P(2,2))}{square root over(P(1,1)+P(2,2))}

Similarly, the Vertical Figure of Merit may be computed by the processor16 as the 2-sigma vertical error from the error covariance matrix

VFOM=2√{square root over (P(3,3))}

Parity Space RAIM

At a step 250, the processor 16 computes at least one protection levelvalue. In doing so, the processor 16 can employ a parity spacetechnique. Recall that the weighted least squares solution may beexpressed as:

$\begin{matrix}{{\Delta \hat{x}} = {\left( {H^{T}{WH}} \right)^{- 1}H^{T}W\; \Delta \overset{\sim}{\rho}}} \\{= {S\; \Delta \overset{\sim}{\rho}}} \\{W = {R^{- 1} = \left( {E\left\lbrack {\delta\rho\delta\rho}^{T} \right\rbrack} \right)^{- 1}}}\end{matrix}$

The correlated measurement set can be transformed by the processor 16into a set that is uncorrelated by factoring, using a matrixfactorization method known in the art, the weighting matrix W:

W=L^(T)L

Where L is the lower triangular square root of W. This results in:

$\begin{matrix}{{\Delta \hat{x}} = {{\left( {H^{T}L^{T}{LH}} \right)^{- 1}H^{T}L^{T}L\; \Delta \overset{\sim}{\rho}} = {\left\lbrack {({LH})^{T}{LH}} \right\rbrack^{- 1}({LH})^{T}L\; \Delta \overset{\sim}{\rho}}}} \\{= {{\left( {{\overset{\_}{H}}^{T}\overset{\_}{H}} \right)^{- 1}{\overset{\_}{H}}^{T}\Delta \overset{\_}{\rho}} = {\overset{\_}{S}\; \Delta \overset{\_}{\rho}}}}\end{matrix}$${{{where}\mspace{14mu} \overset{\_}{H}} = {LH}},{{\Delta \overset{\_}{\rho}} = {L\; \Delta \overset{\sim}{\rho}}},{\overset{\_}{S} = {\left( {{\overset{\_}{H}}^{T}\overset{\_}{H}} \right)^{- 1}\overset{\_}{H}}}$

The covariance of the transformed measurement errors is:

$\begin{matrix}{{E\left\lbrack {\delta \overset{\_}{\rho}\delta {\overset{\_}{\rho}}^{T}} \right\rbrack} = {E\left\lbrack {L\; {\delta\rho\delta\rho}^{T}L^{T}} \right\rbrack}} \\{= {{{LW}^{- 1}L^{T}} = {{L\left( {L^{T}L} \right)}^{- 1}L^{T}}}} \\{= {{L\left( {L^{- 1}L^{- T}} \right)}L^{T}}} \\{= I_{N}}\end{matrix}$

Thus, it is seen that the transformed measurements are uncorrelated,each with unity variance. Multiplying our measurement equation by L theprocessor 16 gets:

LΔ{tilde over (ρ)}=LHΔx+Lδρ

Δ ρ= HΔx+δ ρ

An N×N orthogonal matrix Q can be found by the processor 16, such that:

${Q\overset{\_}{H}} = {{\begin{bmatrix}q_{11} & q_{12} & \ldots & q_{1N} \\q_{21} & q_{22} & \ldots & q_{2N} \\q_{31} & q_{32} & \ldots & q_{3N} \\q_{41} & q_{42} & \ldots & q_{4N} \\q_{51} & q_{52} & \ldots & q_{5N} \\\vdots & \vdots & \vdots & \vdots \\q_{N\; 1} & q_{N\; 2} & q_{N\; 3} & q_{NN}\end{bmatrix}\overset{\_}{H}} = {{\begin{bmatrix}A \\B\end{bmatrix}\overset{\_}{H}} = {\begin{bmatrix}u_{11} & u_{12} & u_{13} & u_{14} \\0 & u_{22} & u_{23} & u_{24} \\0 & 0 & u_{33} & u_{34} \\0 & 0 & 0 & u_{44} \\0 & 0 & 0 & 0 \\\vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & 0\end{bmatrix} = \begin{bmatrix}U \\0\end{bmatrix}}}}$

Thus, if the processor 16 multiplies the measurement equation by Q, theresult is:

${Q\; \Delta \overset{\_}{\rho}} = {{{Q\overset{\_}{H}\Delta \; x} + {Q\; \delta {\overset{\_}{\rho}\begin{bmatrix}A \\B\end{bmatrix}}\Delta \overset{\_}{\rho}}} = {{\begin{bmatrix}U \\0\end{bmatrix}\Delta \; x} + {\begin{bmatrix}A \\B\end{bmatrix}\delta \overset{\_}{\rho}}}}$

Where A is the upper 4×N portion of Q and B is the lower (N−4)×Nportion. The following two equations result:

AΔ ρ= UΔx+Aδ ρ

BΔ ρ=Bδ ρ=p

The first equation can be used to solve for the estimated least squaressolution by setting the pseudo-range error to zero:

AΔ ρ=UΔ{circumflex over (x)}

Δ{circumflex over (x)}=U⁻¹AΔ ρ

This equation is a more efficient than the previously given one since itonly requires inverting an upper-triangular matrix.

The second equation gives us the parity vector p which is non-zero onlywhen there are pseudo-range errors. In the absence of a failure, theparity covariance is:

$\begin{matrix}{{E\left\lbrack {pp}^{T} \right\rbrack} = {{{{BE}\left\lbrack {\delta \overset{\_}{\rho}\delta {\overset{\_}{\rho}}^{T}} \right\rbrack}B^{T}} = {BB}^{T}}} \\{= I_{N - 4}}\end{matrix}$

Thus the parity elements are also uncorrelated zero mean Gaussian randomvariables with unity variance.

In an embodiment, the processor 16 employs a chi-square method using theconcept of pbias. In such an embodiment, the processor 16 uses thesquare of the parity magnitude as the discriminator (test statistic) das follows:

d=p ^(T) p=p ₁ ² +p ₂ ² + . . . +p _(N-4) ²

The discriminator will then have a central chi-square distribution withN−4 degrees of freedom. The processor 16 places a threshold on thisdiscriminator above which a failure is declared. This threshold D iscomputed by the processor 16 from the chi-square probability densityfunction to yield an allowable false alarm probability.

Once the threshold has been set, the question becomes, with a failure ona single satellite that results in the discriminator just reaching thethreshold, how large the resulting position error can be while meetingthe probability of missed detection. A bias failure ε on the kthsatellite will result in a solution error and a parity magnitude errorof:

${\delta \hat{x}} = {{S\; {\delta\rho}} = {{S\begin{bmatrix}0 \\\vdots \\ɛ \\\vdots \\0\end{bmatrix}} = {\begin{bmatrix}\begin{matrix}\begin{matrix}s_{1\; k} \\s_{2\; k}\end{matrix} \\s_{3\; k}\end{matrix} \\s_{4\; k}\end{bmatrix}ɛ}}}$ $\begin{matrix}{{p} = \sqrt{p^{T}p}} \\{= \sqrt{{\delta\rho}^{T}L^{T}B^{T}{BL}\; {\delta\rho}}} \\{= {ɛ\sqrt{\Lambda_{kk}}}} \\{{{where}\mspace{14mu} \Lambda} = {{{\overset{\_}{B}}^{T}\overset{\_}{B}\mspace{14mu} {and}\mspace{14mu} \overset{\_}{B}} = {BL}}}\end{matrix}$

Thus the parity bias is related to the horizontal position error throughthe following slope:

${{slope}(k)} = \frac{\sqrt{s_{1k}^{2} + s_{2k}^{2}}}{\sqrt{\Lambda_{kk}}}$

The satellite with the largest slope is the most difficult to detect.This slope is referred to as Slope_(max). FIG. 3 depicts a noise scatterthat would occur if there was a bias on the most difficult to detectsatellite in the presence of the expected noise. The specific bias thatresults in the percentage of data to the left of the detection thresholdD equal to the missed detection probability is of particular interest.The parity magnitude associated with this bias is called “pbias.”

With the bias present, the discriminator (square of the paritymagnitude) has a non-central chi-square distribution with N−4 degrees offreedom. It can be shown that the non-centrality parameter λ of thechi-square distribution is:

λ=pbias²

Thus, using the non-central chi-square probability density function, theprocessor 16 can determine the value for pbias which meets the requiredprobability of missed detection.

The Horizontal Protection Level (HPL) is then:

HPL=Slope_(max) ·pbias

Gaussian Method by Rotation of the Parity Space

With reference to FIG. 4, through additional orthogonal transformations,an embodiment can rotate the parity space such that the parity error dueto a bias failure on satellite k is entirely along axis 1 of the parityspace.

$\begin{matrix}{p^{k} = {Q^{k}\overset{\_}{B}{\delta\rho}}} \\{= {{\overset{\_}{B}}^{k}{\delta\rho}}} \\{= {\begin{bmatrix}{\overset{\_}{b}}_{1}^{k} & \ldots & {\overset{\_}{b}}_{k}^{k} & \ldots & {\overset{\_}{b}}_{N}^{k} \\X & \vdots & 0 & \vdots & X \\\vdots & \vdots & \vdots & \vdots & \vdots \\X & \vdots & 0 & \vdots & X\end{bmatrix}\begin{bmatrix}0 \\\vdots \\ɛ \\\vdots \\0\end{bmatrix}}} \\{= {\begin{bmatrix}\left( {\overset{\_}{b}}^{k} \right)^{T} \\X \\X \\X\end{bmatrix}\begin{bmatrix}0 \\\vdots \\ɛ \\\vdots \\0\end{bmatrix}}} \\{= {{\overset{\_}{b}}_{k}^{k}ɛ}}\end{matrix}$

Since the bias only shows up on axis 1, the result is a scalar, and thediscriminator, in general, can be determined by the processor 16 usingthe following:

d ^(k) =p ₁ ^(k)=( b ^(k))^(T)δρ

With no failure, only correlated random errors w are present:

d^(k)=( b ^(k))^(T) w=(b^(k))^(T) Lw=(b^(k))^(T) w

where B^(k)=Q^(k)B= B ^(k)L⁻¹, w=Lw, E[ w w ^(T)]=I_(N)

Note: σ_(d) _(k) ²=E[(d^(k))²]=E[(b^(k))^(T) wwb^(k)]=E[(b^(k))^(T)b^(k)]=1

With a bias failure on satellite k plus the correlated random errors won each of the satellites, the discriminator can be determined by theprocessor 16 using the following:

d ^(k) = b _(k) ^(k)ε+(b ^(k))^(T) w

The impact of the noise in addition to failure on the horizontalposition is:

$\begin{matrix}{{\delta \; r^{h}} = {{S^{h}\delta \; \rho} = {{{s_{k}^{h}ɛ} + {S^{h}w}} = {{s_{k}^{h}ɛ} + {{\overset{\_}{S}}^{h}{Lw}}}}}} \\{= {{s_{k}^{h}ɛ} + {{\overset{\_}{S}}^{h}\overset{\_}{w}}}}\end{matrix}$

where S^(h), S ^(h)=1^(st) two rows of S and S, and s_(k) ^(h)=column kof S^(h)

Using the Gaussian probability density function, the threshold D, whichmeets the probability of false alarm, can be determined by the processor16. At detection, the discriminator magnitude is equal to the threshold:

|d ^(k) |=| b ^(k)ε+(b ^(k))^(T) w|=D

Assuming the failure is positive and much larger than the noise:

${{{{\overset{\_}{b}}_{k}^{k}}ɛ} + {\left( b^{k} \right)^{T}\overset{\_}{w}}} = {\left. D\Rightarrow ɛ \right. = {\frac{D}{{\overset{\_}{b}}_{k}^{k}} - \frac{\left( b^{k} \right)^{T}\overset{\_}{w}}{{\overset{\_}{b}}_{k}^{k}}}}$

The resulting horizontal position error can be determined by theprocessor 16 using the following:

$\begin{matrix}{{\delta \; r^{h}} = {{{s_{k}^{h}ɛ} + {{\overset{\_}{S}}^{h}\overset{\_}{w}}} = {{s_{k}^{h}\left( {\frac{D}{{\overset{\_}{b}}_{k}^{k}} - \frac{\left( b^{k} \right)^{T}\overset{\_}{w}}{{\overset{\_}{b}}_{k}^{k}}} \right)} + {{\overset{\_}{S}}^{h}\overset{\_}{w}}}}} \\{= {{s_{k}^{h}\left( {\frac{D}{{\overset{\_}{b}}_{k}^{k}} - {\frac{1}{{\overset{\_}{b}}_{k}^{k}}{\sum\limits_{i = 1}^{N}{b_{i}^{k}{\overset{\_}{w}}_{i}}}}} \right)} + {\sum\limits_{i = 1}^{N}{{\overset{\_}{s}}_{i}^{h}{\overset{\_}{w}}_{i}}}}}\end{matrix}$

The position error magnitude in the direction of the failure is:

$\begin{matrix}{{\delta \; r_{k}^{h}} = {{\frac{s_{k}^{h}}{s_{k}^{h}} \cdot \delta}\; r^{h}}} \\{= {{{s_{k}^{h}}\left( {\frac{D}{{\overset{\_}{b}}_{k}^{k}} - {\frac{1}{{\overset{\_}{b}}_{k}^{k}}{\sum\limits_{i = 1}^{N}{b_{i}^{k}{\overset{\_}{w}}_{i}}}}} \right)} + {\sum\limits_{i = 1}^{N}{\frac{s_{k}^{h} \cdot {\overset{\_}{s}}_{i}^{h}}{s_{k}^{h}}{\overset{\_}{w}}_{i}}}}} \\{= {\underset{\underset{Bias}{}}{\frac{s_{k}^{h}}{{\overset{\_}{b}}_{k}^{k}}D} + \underset{\underset{Noise}{}}{\frac{1}{{\overset{\_}{b}}_{k}^{k}}{\sum\limits_{i = 1}^{N}{\left( {{\frac{s_{k}^{h} \cdot {\overset{\_}{s}}_{i}^{h}}{s_{k}^{h}}{{\overset{\_}{b}}_{k}^{k}}} - {{s_{k}^{h}}b_{i}^{k}}} \right){\overset{\_}{w}}_{i}}}}}}\end{matrix}$

Since the random pseudo-range error is uncorrelated with unity variance,the variance of the noise term about the mean can be determined by theprocessor 16 using the following:

$\sigma_{\delta \; r_{k}^{h}} = \sqrt{\frac{1}{{\overset{\_}{b}}_{k}^{k}}{\sum\limits_{i = 1}^{N}\left( {{\frac{s_{k}^{h} \cdot {\overset{\_}{s}}_{i}^{h}}{s_{k}^{h}}{{\overset{\_}{b}}_{k}^{k}}} - {{s_{k}^{h}}b_{i}^{k}}} \right)^{2}}}$

The Horizontal Protection Level for satellite k can be determined by theprocessor 16 using the following:

$\begin{matrix}{{HPL}_{k} = {{\frac{s_{k}^{h}}{s_{k}^{h}}D} + {K_{md}\sigma_{\delta \; r_{k}^{h}}}}} \\{= {{\frac{s_{k}^{h}}{s_{k}^{h}}K_{fa}} + {K_{md}\sigma_{\delta \; r_{k}^{h}}}}}\end{matrix}$

where K_(md) and K_(fa) are sigma multipliers set to meet theprobabilities of false alarm and missed detection

This process is repeated for all N satellites and the total HPL can bedetermined by the processor 16 using the following:

HPL=max(HPL_(k)), k=1,N

Ionospheric Error Model Calculations

Determination of Ionospheric Grid Points and Pierce Point Coordinates

For a Kalman filter approach, and in order to utilize (17), theprocessor 16 may first determine the coordinates of each gridpoint andthe coordinates of the satellite's ionospheric pierce point. Then, usingthose two sets of coordinates, the great circle distance between thepierce point and the grid point can be calculated by the processor 16.For either a Kalman filter or snapshot RAIM approach, knowing thecoordinates of a point i (e.g., the system illustrated in FIG. 1, or“user”) and the distance and azimuth from point i to a point j (e.g. thegridpoint), the coordinates of point j can be determined by theprocessor 16 as follows:

$\begin{matrix}{\lambda_{j} = {\tan^{- 1}\left( \frac{{\sin \; \lambda_{i}\cos \; \psi_{ij}} + {\cos \; \lambda_{i}\sin \; \psi_{ij}\cos \; A_{ij}}}{\sqrt{\left( {{\cos \; \lambda_{i}\cos \; \psi_{ij}} - {\sin \; \lambda_{i}\sin \; \psi_{ij}\cos \; A_{ij}}} \right)^{2} + {\sin^{2}\psi_{ij}\sin^{2}A_{ij}}}} \right)}} & \left( {A{.1}} \right) \\{\mspace{79mu} {{{\Lambda_{j} = {\Lambda_{i} + {\tan^{- 1}\left( \frac{\sin \; \psi_{ij}\sin \; A_{ij}}{{\cos \; \lambda_{i}\cos \; \psi_{ij}} - {\sin \; \lambda_{i}\sin \; \psi_{ij}\cos \; A_{ij}}} \right)}}}\mspace{79mu} {\lambda_{i} = {{Geodetic}\mspace{14mu} {latitude}\mspace{14mu} {of}\mspace{14mu} {point}\mspace{14mu} i}}\mspace{79mu} {\lambda_{j} = {{Geodetic}\mspace{14mu} {latitude}\mspace{14mu} {of}\mspace{14mu} {point}\mspace{14mu} j}}\mspace{79mu} {\Lambda_{i} = {{Geodetic}\mspace{14mu} {longitude}\mspace{14mu} {of}\mspace{14mu} {point}\mspace{14mu} i}}\mspace{79mu} {\Lambda_{j} = {{Geodetic}\mspace{14mu} {longitude}\mspace{14mu} {of}\mspace{14mu} {point}\mspace{14mu} j}}\mspace{79mu} {A_{ij} = {{Azimuth}\mspace{14mu} {angle}\mspace{14mu} ({bearing})\mspace{14mu} {from}\mspace{14mu} {point}\mspace{14mu} i\mspace{14mu} {to}\mspace{14mu} {point}\mspace{14mu} j}}{\psi_{ij} = {{{Angular}\mspace{14mu} {distance}\mspace{14mu} \left( {{earth}^{\prime}s\mspace{14mu} {central}\mspace{14mu} {angle}} \right)\mspace{20mu} {from}\mspace{14mu} {point}\mspace{14mu} i\mspace{14mu} {to}\mspace{14mu} {point}\mspace{14mu} j}\; \mspace{104mu} = \frac{d_{ij}}{R_{c} + h_{l}}}}}\mspace{79mu} {d_{ij} = {{Great}\mspace{14mu} {circle}\mspace{14mu} {distance}\mspace{14mu} {from}\mspace{14mu} {point}\mspace{14mu} i\mspace{14mu} {to}\mspace{14mu} {point}\mspace{14mu} j}}\mspace{79mu} {R_{e} = {{{Radius}\mspace{14mu} {of}\mspace{14mu} {the}\mspace{14mu} {earth}} = {6378\mspace{14mu} {km}}}}\mspace{79mu} {h_{l} = {{{Height}\mspace{14mu} {of}\mspace{14mu} {the}\mspace{14mu} {ionosphere}\mspace{14mu} {thin}\mspace{14mu} {shell}\mspace{14mu} {model}} = {350\mspace{14mu} {km}}}}}} & \left( {A{.2}} \right)\end{matrix}$

The coordinates of the ionospheric pierce point of the satellite canalso be calculated using (A.1) and (A.2). In this case, ψ_(ij)represents the central angle from the user location to the pierce pointand may be calculated by the processor 16 as follows:

$\begin{matrix}{\psi_{ij} = {\frac{\pi}{2} - E - {\sin^{- 1}\left( {\frac{R_{e}}{R_{e} + h_{l}}\cos \; E} \right)}}} & \left( {A{.3}} \right)\end{matrix}$

where E is the elevation angle of the satellite from the user locationwith respect to the local tangent plane.

Computing Elevation and Azimuth Angles of Satellite

The elevation angle E of a satellite is defined as the angle theline-of-sight vector makes with the user's local tangent (horizontal)plane. The azimuth angle A of the satellite is the angle of theline-of-sight vector with respect to true north as measured in thehorizontal plane. Thus, we have the following

$\begin{matrix}{E = {{A{TAN}}\; 2\left( {{- u_{LOS\_ z}},\sqrt{u_{LOS\_ x}^{2} + u_{LOS\_ y}^{2}}} \right)}} & \left( {A{.4}} \right) \\{A_{0} = {{{A{TAN}}\; 2\left( {u_{LOS\_ y},u_{LOS\_ x}} \right)} + \alpha}} & \left( {A{.5}} \right) \\{A = \left\{ \begin{matrix}{A_{0},} & {{- \pi} \leq A_{0} < \pi} \\{{A_{0} - {2\; \pi}},} & {A_{0} \geq \pi} \\{{A_{0} + {2\; \pi}},} & {A_{0} < \pi}\end{matrix} \right.} & \left( {A{.6}} \right)\end{matrix}$

where

u_(LOS) _(—) _(x),u_(LOS) _(—) _(y),u_(LOS) _(—) _(z)=x, y, and zcomponents of the line-of-sight vector u_(LOS) ^(L)

α=wander angle (angle in azimuth from north to the x local-level frameaxis)

Note that the azimuth angle is adjusted by ±2π so that the result isbetween −π and +π.

Determination of Great Circle Distance

The great circle distance along the ionospheric thin shell model from apoint i (e.g. satellite pierce point) to another point j (e.g. gridpoint) may be calculated by the processor 16 as follows:

$\begin{matrix}{d_{ij} = {\left( {R_{e} + h_{l}} \right){\tan^{- 1}\left( \frac{\sqrt{{\cos^{2}\lambda_{j}\sin^{2}{\Delta\Lambda}_{ij}} + \left( {{\cos \; \lambda_{i}\sin \; \lambda_{j}} - {\sin \; \lambda_{i}\cos \; {\Delta\Lambda}_{ij}}} \right)^{2}}}{{\cos \; \lambda_{i}\cos \; \lambda_{j}\cos \; {\Delta\Lambda}_{ij}} + {\sin \; \lambda_{i}\sin \; \lambda_{j}}} \right)}}} & \left( {A{.7}} \right)\end{matrix}$

where

ΔΛ_(ij)=Λ_(j)−Λ_(i)

Ionospheric Variance Model

The algorithm that may be executed by the processor 16 for calculationof the ionospheric model error variance may be from ICD-GPS-200C andDO-229D J.2.3. Note that the symbols in this section are unique to thissection.

Using the satellite's elevation angle E, form the earth's central anglebetween the user position and the earth projections of ionosphericpierce point ψ_(pp) using equation (A.3).

Next, using the satellite's elevation angle E, the azimuth angle A, theearth's central angle ψ_(pp) and the user geodetic latitude λ_(u) andlongitude Λ_(u), determine the pierce point geodetic latitude φ_(pp) andlongitude λ_(pp) using equations (A.1) and (A.2).

Form the absolute value of the geomagnetic latitude of the ionosphericpierce point.

|λ_(m)|=|λ_(pp)+0.064π cos(Λ_(pp)−1.617π)|radians  (A.8)

Form an estimate of the vertical delay error based on geomagneticlatitude

$\begin{matrix}{\tau_{vert} = \left\{ \begin{matrix}{{9\mspace{14mu} {meters}},} & {{\lambda_{m}} \leq {20\mspace{14mu} {degrees}}} \\{{4.5\mspace{14mu} {meters}},} & {{22.5{^\circ}} < {\lambda_{m}} \leq {55{^\circ}}} \\{{6\mspace{14mu} {meters}},} & {{\lambda_{m}} > {55{^\circ}}}\end{matrix} \right.} & \left( {A{.9}} \right)\end{matrix}$

Using the elevation angle E, calculate the square of the obliquityfactor.

$\begin{matrix}{F_{pp}^{2} = \frac{1}{1 - \left( \frac{R_{e}{\cos (E)}}{R_{e} + h_{l}} \right)^{2}}} & \left( {A{.10}} \right)\end{matrix}$

Form the modeled estimated variance of the ionospheric delay.

σ_(model) ²=F_(pp) ²τ_(vert) ²  (A.11)

Form the estimated variance using the compensation that is applied ifavailable from the receiver. (If not, assume zero).

$\begin{matrix}{\sigma_{comp}^{2} = \left( \frac{{CT}_{IONO}}{5} \right)^{2}} & \left( {A{.12}} \right)\end{matrix}$

Form the estimated variance of the ionospheric delay.

σ_(iono) ²=max(σ_(model) ²,σ_(comp) ²)  (A.13)

While the preferred embodiment of the invention has been illustrated anddescribed, as noted above, many changes can be made without departingfrom the spirit and scope of the invention. Accordingly, the scope ofthe invention is not limited by the disclosure of the preferredembodiment. Instead, the invention should be determined entirely byreference to the claims that follow.

1. A navigation system for a vehicle having a receiver operable toreceive a plurality of signals from a plurality of transmitters, thenavigation system comprising: a processor; and a memory device havingstored thereon machine-readable instructions that, when executed by theprocessor, enable the processor to: determine a set of error estimatescorresponding to pseudo-range measurements derived from the plurality ofsignals, determine an error covariance matrix for a main navigationsolution using ionospheric-delay data, and using a parity spacetechnique, determine at least one protection level value based on theerror covariance matrix.
 2. The system of claim 1 wherein determiningthe error covariance matrix includes determining a spatially correlatedionospheric error associated with each of the transmitters.
 3. Thesystem of claim 2 wherein determining the error covariance matrixincludes defining a plurality of grid points on a thin shell model ofthe ionosphere.
 4. The system of claim 3 wherein the defined grid pointsare equally spaced in azimuth at a great circle distance from thesystem.
 5. The system of claim 4 wherein the great circle distance is1500 km.
 6. The system of claim 3 wherein determining the errorcovariance matrix includes defining a vector of normalized ionosphericdelays at the grid points.
 7. The system of claim 1 wherein the at leastone protection level value comprises a horizontal protection levelvalue.
 8. The system of claim 1 wherein the at least one protectionlevel value comprises a vertical protection level value.
 9. The systemof claim 1 wherein the at least one protection level value comprises ahorizontal exclusion level value.
 10. A computer-readable medium havingcomputer-executable instructions for performing steps comprising:determining a set of error estimates corresponding to pseudo-rangemeasurements derived from the plurality of signals; determining an errorcovariance matrix for a main navigation solution using ionospheric-delaydata; using a parity space technique, determining at least oneprotection level value based on the error covariance matrix; anddisplaying the at least one protection level value.
 11. The medium ofclaim 10 wherein determining the error covariance matrix includesdetermining a spatially correlated ionospheric error associated witheach of the transmitters.
 12. The medium of claim 11 wherein determiningthe error covariance matrix includes defining a plurality of grid pointson a thin shell model of the ionosphere.
 13. The medium of claim 12wherein the defined grid points are equally spaced in azimuth at a greatcircle distance from the system.
 14. The medium of claim 13 wherein thegreat circle distance is 1500 km.
 15. The medium of claim 12 whereindetermining the error covariance matrix includes defining a vector ofnormalized ionospheric delays at the grid points.
 16. The medium ofclaim 10 wherein the at least one protection level value comprises ahorizontal protection level value.
 17. The medium of claim 10 whereinthe at least one protection level value comprises a vertical protectionlevel value.
 18. The medium of claim 10 wherein the at least oneprotection level value comprises a horizontal exclusion level value. 19.A method, comprising the steps of: accessing from a first computer thecomputer-executable instructions of claim 10; and providing theinstructions to a second computer over a communications medium.